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Abstract — The performance of the generalized belief propa- 
gation algorithm to compute the noiseless capacity and mu- 
tual information rates of finite-size two-dimensional and three- 
dimensional run-length limited constraints is investigated. In 
both cases, the problem is reduced to estimating the partition 
function of graphical models with cycles. The partition function 
is then estimated using the region-based free energy approxi- 
mation technique. For each constraint, a method is proposed 
to choose the basic regions and to construct the region graph 
which provides the graphical framework to run the generalized 
belief propagation algorithm. Simulation results for the noiseless 
capacity of different constraints as a function of the size of 
the channel are reported. In the cases that tight lower and 
upper bounds on the Shannon capacity exist, convergence to the 
Shannon capacity is discussed. For noisy constrained channels, 
simulation results are reported for mutual information rates as 
a function of signal-to-noise ratio. 

Index Terms — Generalized belief propagation algorithm, run- 
length limited constraints, partition function, factor graphs, 
region graphs, noiseless capacity, Shannon capacity, mutual 
information rate. 

I. Introduction 

Run-length limited (RLL) constraints are widely used in 
magnetic and optical recording systems. Such constraints 
reduce the effect of inter-symbol interference and help in 
timing control. In track-oriented storage systems constraints 
are defined in one dimension. 

We say a binary one-dimensional (1-D) sequence satisfies 
the (d, k)-RLL constraint if the runs of O's have length at most 
k and the runs of O's between successive l's have length at 
least d. We suppose that < d < k < oo. 

The Shannon capacity of a 1-D (d, fc)-RLL constraint is 
defined as 

C W0 a Hm log 2 Z(m) 
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where Z(m) denotes the number of binary 1-D sequences of 
length m that satisfy the (d, fc)-RLL constraint, see Qj, j2). 

With the rise in demand for larger storage in smaller 
size and with recent developments in page-oriented storage 
systems, such as holographic data storage, two-dimensional 
(2-D) constraints have become more of interest |3|. In these 
systems, data is organized on a surface and constraints are 
defined in two dimensions. 

A 2-D binary array satisfies the {d\, k\, d2, fc2)-RLL con- 
straint if it satisfies a (<f 1; fcjJ-RLL constraint horizontally and 
a (<i2, &2)-RLL constraint vertically. If a 2-D binary array 
satisfies a 1-D (d, fe)-RLL constraint both horizontally and 
vertically, we simply say that it satisfies a 2-D (d, fc)-RLL 
constraint. 

Example: 2-D (2, oo)-RLL constraint: 

The 2-D (2, oo)-RLL constraint is satisfied in the following 
2-D binary array segment. In words, in every row and every 
column of the array there are at least two O's between succes- 
sive l's; but the runs of O's can be of any length (however, 
l's can be diagonally adjacent). 

. . . 0100100001001000100000100010 . . . 
. . . 1000010000100010000100000100 . . . 
. . . 0001000010000001000000010001 . . . 
. . . 0100100100010000001000100000 . . . 

The Shannon capacity of a 2-D (di, ki, di, fc2)-RLL con- 
straint is defined as 

c £, kM2) . Um log a Z( ro ,n) 
m,n->oo run 

where Z(m,n) denotes the number of 2-D binary arrays of 
size m x n that satisfy the {d\, fei, cfej fc2)-RLL constraint. 

Similarly, the Shannon capacity can be defined for higher 
dimensional constrained channels. For example, the Shannon 
capacity in three dimensions (J^' kl ' d2 ' k2 < d3 ' k3 ) depends on 
Z(m, n, q), the number of three-dimensional (3-D) binary 
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arrays of size m x n x q that satisfy a {d\, ki,d,2, fe, ^3, fe)- 
RLL constraint. 

The noiseless capacity of a constrained channel is an impor- 
tant quantity that provides an upper bound to the information 
rate of any encoder that maps arbitrary binary input into binary 
data that satisfies a given constraint. There are a number of 
techniques to compute the 1-D Shannon capacity (for example 
combinatorial or algebraic approaches) fT). In contrast to the 
1-D capacity, except for a few cases, exact values of two 
and higher dimensional (positive) Shannon capacities are not 
known, see (4|-[9|. 

For noisy 1-D constrained channels, simulation-based tech- 
niques proposed in (16) , (T7| can be used to compute mutual 
information rates. However, computing mutual information 
rates of noisy 2-D RLL constraints has been an unsolved 
problem. 

In this paper, the goal is to apply the generalized belief 
propagation (GBP) algorithm (10) for the above-mentioned 
problems, namely, to compute an estimate of the capacity of 
noiseless 2-D and 3-D RLL constrained channels and mutual 
information rates of noisy 2-D constrained channels. For both 
problems GBP turns out to yield very good approximate 
results. 

Preliminary versions of the material of this paper have 
appeared in ||TTJ and 1 12 1. In 1 1 1 1, we applied GBP to compute 
the noiseless capacity of 2-D and 3-D RLL constrained chan- 
nels. In (12) , GBP was applied to compute mutual information 
rates of a 2-D (1, oo)-RLL constrained channel with relatively 
small size and only at high signal-to-noise ratio (SNR). In this 
paper, we show that both problems reduce to estimating the 
partition function of graphical models with cycles. We then 
apply GBP to both problems and consider new constraints 
and larger sizes of grid. 

Our main motivations for this research were the successful 
application of GBP for information rates of 2-D finite-state 
channels with memory in [13], Kikuchi approximation for 
decoding of LDPC codes and partial-response channels in 1 14|, 
and tree-based Gibbs sampling for the noiseless capacity and 
information rates of 2-D constrained channels in (12), (15). 

The outline of the paper is as follows. In Section [II] we 
consider the problem of computing the partition function and 
discuss how this problem is related to computing the noiseless 
capacity and information rates of constrained channels. Region 
graphs, GBP, and region-based free energy are outlined in 



capacity of noiseless 2-D and 3-D RLL constraints are reported 
in Section IV-A In Section [V] we apply GBP to compute 
mutual information rates of noisy 2-D RLL constraints and 
report numerical experiments for mutual information rates in 
Section IV-Al 

II. Problem Set-up 

Consider a 2-D channel of size N = m X m with a set 
of X = {X±,X2, ■ ■ ■ ,Xn} random variables. Let Xi denote 
a realization of and let x denote {xi,X2, ■ ■ ■ ,xn}- We 
assume that each takes values in a finite set Xi. Also let 
X be the Cartesian product X = X\ x X2 x . . . x X^. 

In constrained channels, not all sequences of symbols from 
the channel alphabet X are admissible. Let S% C X be the set 
of admissible input sequences. We define the indicator function 

f l, xeS x 
0, x^S x 



(3) 



The partition function Z is defined as 



(4) 



With the above definitions, Z = \S%\ is the number of 
sequences that satisfy a given constraint. Therefore, computing 
the capacity of constrained channels as expressed in |2]), is 
closely related to computing the partition function as in Q. 

Also note that with the above definitions 

*M - ^ (5) 

is a probability mass function on X. 

For a noisy 2-D channel, let X be the input and 
Y = {Yi, Y2, ■ ■ ■ , ijv} be the output of the channel. The 
mutual information rate is 



i/(X;Y) = i(ff(Y)-H(Y|X)). 



(6) 



Let us suppose that H(Y\X) is analytically available. In 
this case, the problem of estimating the mutual information 
rate reduces to estimating the entropy of the channel output, 
which is 

ff(Y) = -E[logp(Y)]. (7) 

As in [16], we can approximate the expectation in (|7]i by 
drawing L samples y^^y*- 2 ', . . . ,y^ L ' according to p(y) and 
use the empirical average as 



Section III Section IV discusses the capacity of noiseless 2-D 
constraints. Numerical values and simulation results for the 



H(Y) 
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Therefore, the problem of estimating the mutual information 
rate reduces to computing p{y^) for I = 1, 2, . . . , L. 
We will compute p(y^') based on 



>(y w ) = 5>(*)p(y w W. 



(9) 



which for a fixed y^ has also the form Q and therefore 
requires the computation of a partition function. 

RLL constraints impose restrictions on the values of vari- 
ables that can be verified locally. For example, in a 2-D (1. oo)- 
RLL constraint no two (horizontally or vertically) adjacent 
variables can both have the value 1. The indicator function of 
this constraint factors into a product of kernels of the form 



0, if Xi = Xj = 1 

1 , else, 



(10) 



with one such kernel for each adjacent pair (xi,Xj). 

The factorization with kernels as in ( fTO) can be represented 
with a graphical model. In this paper, we focus on graphical 
models defined in terms of Forney factor graphs. Fig. [T] shows 
the Forney factor graph of a 2-D (1, oo)-RLL constraint where 
the boxes labeled "=" are equality constraints fl9) . (Fig. [T] 
may also be viewed as a factor graph as in fl8) where the 
boxes labeled "=" are the variable nodes). 

In general, we suppose that the indicator function /(x) of 
an RLL constraint factors into a product of non-negative local 
kernels each having some subset of x as arguments; i.e. 
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where x a is a subset of x and each kernel f a (x a ) has elements 
of x a as arguments. 

In this case, the partition function in Q can be written as 



(12) 



x£X a 



If the factorization in ( fTT) yields a cycle-free factor graph 
(with not too many states), the sum in (12) , or equivalently the 
sum in Q, can be computed efficiently by the sum-product 
message passing algorithm |18]. However, for the examples 
we study in this paper, like the Forney factor graph of a 2-D 
(l,oo)-RLL constraint in Fig. [I] factor graphs contain (many 
short) cycles. In such cases computing Z requires a sum with 
an exponential number of terms and therefore we are interested 
in applying approximate methods. 

Due to the presence of many short cycles in the factor graph 
representation of 2-D and 3-D RLL constraints, loopy belief 
propagation often fails to converge. As a result, we apply GBP 
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Fig. 1. Forney factor graph for a 2-D (1, oo)-RLL constraint. 

to estimate Z, which then leads to estimating the noiseless 
capacity and mutual information rates of RLL constraints. 

III. GBP AND THE REGION GRAPH METHOD 

In statistical physics, Z defined in Q is known as the 
partition function and the Helmholtz free energy is defined 



MZ). 



(13) 



The partition function and the Helmholtz free energy are 
important quantities in statistical physics since they carry 
information about all thermodynamic properties of a system. 

A number of techniques have been developed in statistical 
physics to approximate the free energy. The method that 
we apply in this paper is known as the region-based free 
energy approximation, in particular we use the cluster variation 
method to select a valid set of regions and counting numbers, 
see fl0| and (20) for more details. 

We start by introducing the region graph representation of 
our problem. Such a region graph will provide a graphical 
framework for GBP algorithm. For each RLL constraint, the 
size of the basic region is chosen based on the constraint 
parameters. For a 2-D (d%, ki, fc2)-RLL constraint with 
finite k\ and k 2 , the width and the height of the basic region 
is chosen as 

W R = h + 1 
H R = k 2 + 1, 

and for the infinite case, the size is chosen as 

W R = di + 1 

Hr = d 2 + l. 
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Fig. 2. Basic region of size 2x2 for a 2-D (1, oo)-RLL constraint. 
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Fig. 3. The region graph for Forney factor graph in Fig. [2] 




Such a choice for the basic regions seems plausible since 
the validity of a given array can be determined by verifying 
the constraints in each region and sliding the basic regions 
along the rows and along the columns of the array. For a 
2-D (l,oo)-RLL constraint, Fig. [2] shows the basic regions 
and Fig. [3] shows the region graph and the counting numbers 
associated with each region. 

After forming the region graph using the cluster variation 
method, we perform GBP on this graph by sending messages 
between the regions while performing exact computations 
inside each region. 

We will need the region-based free energy to estimate the 
number of arrays that satisfy a given constraint. Therefore, 
we operate GBP on the corresponding region graph until 
convergence and use the obtained region beliefs {6r(x^)} 
to compute the region-based free energy Fjj (as an estimate 
of Fu)- The region-based free energy Fjj can then be used to 
estimate the partition function Z using fjj) . We compute Fjj 
as 



F H = mm F R ({b R (xi R )}) 
{bat 

C*X>*(x*)(to&*(Xfl) ~ ln Il /a( X a))( 14 ) 



Ren x R 



aeA T 



Here 1Z denotes the set of all regions, cr is the counting 
number, stands for the set of variables in region R, and 
Ar is the set of factors in region R. See Fig. [3] 

IV. Capacity of Noiseless 2-D RLL Constraints 

For a 2-D RLL constrained channel of width m and of size 
N = to x to, we run GBP on the corresponding region graph 
to compute Fh and an estimate of Z. We can then compute 

log 2 Z(m, to) 



where Z(m, to) denotes the number of 2-D binary arrays of 
size to x to that satisfy the constraint. 

In our numerical experiments in Section |IV-A| for different 
RLL constraints we show convergence of C(m, m) to the 
Shannon capacity as to increases. 

For example, let us consider a 2-D (l,oo)-RLL constraint 
with corresponding Forney factor graph in Fig. [T] For this 
constraint, we chose basic regions with size 2 x 2 in a sliding 
window manner over the factor graph, see Fig. [2] Starting from 
such basic regions, we applied the cluster variation method on 
the factor graph in Fig. [2] to obtain the corresponding region 
graph depicted in Fig. [3] The counting numbers {cr} are 
shown next to each region. 



A. Numerical Experiments 

Here we present the numerical results of applying GBP to 
estimate the finite-sized noiseless capacity of RLL constraints. 

Tight lower and upper bounds were given for the Shannon 
capacity of a 2-D (1, 00) -RLL constraint in Q. The bounds 
were further improved in J22") and f23) , now known to nine 
decimal digits. 



0.5878911617... < C^ oo) < 0.5878911618. 



(16) 



C(to, to) 



TO X TO 



(15) 



For this constraint, Fig. [4] shows C(to, to) defined in ( fT5| ) 
versus the channel width m over the interval [2,300]. The 
estimation was performed using the parent-to-child and two- 
way GBP algorithms. The two algorithms give almost identical 
results. The horizontal line in Fig. |4] shows the Shannon 
capacity for this channel in ([16}. For a channel of width 300, 
the estimated noiseless capacity is about 0.5884. 
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Fig. 4. Estimated capacity (in bits per symbol) vs. channel width m for Fig. 6. Estimated capacity (in bits per symbol) vs. channel width m for a 
a 2-D (1, oo)-RLL constraint. The horizontal dotted line shows the Shannon 2-D (1, oo, 2, 4)-RLL and (1, oo, 2, 3)-RLL constraints, 
capacity for this channel as in 
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Fig. 5. Estimated capacities (in bits per symbol) vs. channel width m for a Fig. 7. Estimated capacity (in bits per symbol) vs. channel width m for a 



class of 2-D (1, oo, d, oo)-RLL constraints with d = (1, 2, 3, 4). 



2-D (2, oo)-RLL constraint. 



Shown in Fig.|5]are plots of C(m, m) for 2-D (1, oo, d, oo)- 
RLL constraints with d = (1,2,3,4) from top to bot- 
tom, versus the channel width m over the interval [2,200]. 
Fig. [6] shows the plots of C(m,m) for 2-D (1, oo, 2,4)- 
RLL and (1, oo, 2, 3) -RLL constraints versus m over the 
interval [4,200]. From our simulation results, for a chan- 
nel of width 200 the estimated noiseless capacities for 2-D 
(1, oo, d, oo)-RLL constraints with d — (2,3,4) are about 
(0.4994,0.4346,0.3864) and the estimated noiseless capaci- 
ties for 2-D (1, oo, 2, 4)-RLL and (1, oo, 2, 3)-RLL are about 
(0.3106, 0.2109). To the best of our knowledge, no theoretical 
upper or lower bounds exist for these constraints. All plots 
are obtained using the parent-to-child algorithm. Note that 
2-D (1, oo, 1, oo)-RLL plot in Fig. [5] is the same as the plot 
in Fig. g] 

Also shown in Fig. [7] is the plot of C(m, m) for a 2-D 
(2, oo)-RLL constraint versus m over the interval [3, 400]. For 




6 8 10 
Channel width > 



Fig. 8. Estimated capacity (in bits per symbol) vs. channel width m for 
a 3-D (l,oo)-RLL constraint. The horizontal dotted lines show upper and 
lower bounds on the Shannon capacity for this channel as in (18). 



a channel of width 400, the estimated noiseless capacity is 
about 0.4462. Best known lower and upper bounds for the 



6 



Shannon capacity of a 2-D (2, oo)-RLL constraint are given 
in 1 8 1 and |9| respectively, as 



0.4453 < C^ oo) < 0.4457 



(17) 



Our proposed method can be generalized to compute the 
noiseless capacity of 3-D and higher dimensional RLL con- 
straints. For a 3-D (l,oo)-RLL constraint the following lower 
and upper bounds were introduced in |23| 

(18) 



0.5225017418... < C^ D °° ] < 0.5268808478. 



Fig. [8] shows the noiseless capacity estimates of a 3- 
D (l,oo)-RLL constraint, obtained using the parent-to-child 
algorithm, versus the channel width to. The horizontal dotted 
lines show the upper and lower bounds for the Shannon 
capacity. For a channel of width m = 40 the GBP estimated 
capacity is about 0.5267 which falls within these bounds. 

Simulation results and numerical values for the noiseless 
capacity of many other 2-D RLL constraints are reported 
in ||2TJ. 

B. Bounds for the Shannon Capacity 

For any finite to, it is possible to compute lower and upper 
bounds on the Shannon (infinite-size) capacity using C(m, to) 
the capacity of a 2-D RLL constrained channel of width to. 

For example, consider a 2-D (l,oo)-RLL constraint with 
local kernels as in dTOt . From tiling the whole plane with 
to x to squares, it is clear that C(m, to) is an upper bound for 
the Shannon capacity . On the other hand, by tiling the 

plane with m x to squares separated by all-zero guard rows 
id all- 

(l,oo) 



and all-zero guard columns, we obtain (^p±) 2 C(m, to) < 



^223 ■ 

From Fig. [4] the estimated capacity at to = 300 is about 
C(300, 300) = 0.5884, we thus obtain the following lower 
and upper bounds for the Shannon capacity 

0.5844 < C^ oo) < 0.5884. 

Note that although GBP performs remarkably well for 2- 
D constrained channels, it is an approximate algorithm which 
yields approximations to the lower and upper bounds to the 
Shannon capacity. However in order to achieve a desired 
precision, the bounds could provide a criterion for choosing 
the value of to. 

V. Information Rates of Noisy 2-D RLL 
Constraints 

As explained in Section|lIJ the problem of computing mutual 
information rates reduces to computing the output probability. 




Fig. 9. Extension of Fig. [ljto a Forney factor graph of p(x, y) with p(y|x) 
as in 122). 



Therefore, the remaining tasks are 

1) Drawing input samples x^^x' 2 ), . . . ,x^ L ' from Sx ac- 
cording to p(x) and therefrom creating output samples 
y (1) ,y (2) , . . . ,y (L) . 

2) Computing p{y lyt> ) for each I = 1, 2, . . . , L. 

We will compute p(y^) based on 

p(y (f) ) = E?»Wp(y w W, d9) 

where p(x) is a probability mass function on Sx- 

Let us assume uniform distribution over the admissible 
channel input configurations. Therefore we have 



(20) 
(21) 



we also assume the channel is memoryless and p(y|x) 
factors as 



A? 



p(yl x ) = Y[p(yi 



(22) 



For such a noisy 2-D constrained channel, the corresponding 
Forney factor graph, as an extension of Fig. [T] is shown in 
Fig- E 



Using (21 1 and (22 1, we can rewrite (19i as 



JV 

p(yW) = 2 -NC 2D £ l[ p( yW\ Xi)) (23 ) 

xes x i= i 

= 2- NC2D Z(y^), (24) 

where Z(y^) has the same form as the sum in ( fl2| ). 

The input samples x' 1 '^ 2 ', . . . ,x( L ) are generated as fol- 
lows. We run GBP on Fig. [3] until convergence to compute 
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the region beliefs {6r(x.r)} at each region R. The region be- 
liefs are GBP approximations to the corresponding marginals 
{pr( x r}}- In our numerical experiments, each sample X" is 
then generated piecewise sequentially according to the beliefs 
&h(Xr) in basic regions. For example, in the region graph 
of Fig. [3] after computing &r(:ei, x%, Xi,x*>), sample x\ is 
drawn according to bji(xi), sample x 2 is drawn according to 
bR{x2\x\), etc. The input samples x^\x^ 2 \ . . . ,x^ are then 
used to create output 



A 2 ) 



,y( L ) using (22 1 



The beliefs are directly proportional to the factor nodes 
involved in each region, which guarantees that the samples 
are drawn from 6>x- Moreover, since beliefs are good approx- 
imations to the marginal probabilities, one expects that the 
samples are drawn from a distribution close to p(x), see flO) . 

In order to compute Z(y^), as in Section III we start from 
the factor graph in Fig.|9]to build the region graph representing 
the problem and run GBP on this region graph. Finally, the 
estimated p{y^),p(y^), ■ ■ ■ ,p(y^) are used to compute an 
estimate of H (Y) as in 

A. Numerical Experiments 

In our numerical experiments we consider (1, oo)-RLL and 
(2, oo)-RLL constrained channels with size N = 30 x 30 and 
input alphabet X = {-1, +1}^. 

Noise is assumed to be i.i.d. zero mean Gaussian with 
variance a 2 and independent of the input. We thus have 



i/(Y|X) = ^log(27rea 2 ) 



and p(yjx) in ( 22 1 has kernels of the form 
1 



P(Vi\xi) = 



SNR is defined as 



: exp 



SNR = 10 log 



10 a 2 



(25) 



(26) 



(27) 



Shown in Fig. [TUJis the estimated information rate vs. SNR 
over the interval [—10,10] dB for a noisy 2-D (l,oo)-RLL 
constraint. The horizontal dotted line shows the estimated 
noiseless capacity which can be read from Fig. |4] and is about 
0.5943 for this size of channel. 



Illustrated in Fig. 11 is the estimated information rate 



vs. SNR over the interval [—10,10] dB for a noisy 2-D 
(2, oo)-RLL channel. The horizontal dotted line shows the 
estimated noiseless capacity which can be read from Fig. [7] 
and is about 0.4552 for this size of channel. 

The simulation results were obtained by averaging over 
L = 1000 realizations of the channel output. 




Fig. 10. Estimated information rate (in bits per symbol) vs. SNR (in dB) for 
a 30 X 30 channel with a (1, oo)-RLL constraint and additive white Gaussian 
noise. 

0.5 1 1 1 , 1 1 1 1 1 1 1 




Fig. 11. Estimated information rate (in bits per symbol) vs. SNR (in dB) for 
a 30 X 30 channel with a (2, oo)-RLL constraint and additive white Gaussian 
noise. 



Simulation results and numerical values for mutual infor- 
mation rates of many other 2-D RLL constraints are reported 
in (3T). 

VI. Concluding remarks 

We proposed a GBP-based method to estimate the noiseless 
capacity and mutual information rates of RLL constraints in 
two and three dimensions. For noiseless RLL constraints, the 
method was applied to estimate the finite-size capacity of 
different constraints and to show convergence to the Shannon 
capacity as the size of the channel increases. In particular, 
the proposed method can be used to estimate the noiseless 
capacity of RLL constraints in the cases that the capacity is 
not known to a useful accuracy. The method was also applied 
to estimate mutual information rates of noisy RLL constraints 
with additive white Gaussian noise and with a uniform distri- 
bution over the admissible input configurations. Our simulation 
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results show mutual information rates of different constraints 
as a function of SNR. 
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